Force heterogeneities in particle assemblies: From order to disorder 
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The effect of increasing structural disorder on the distribution of contact forces P(f), inside 
three dimensional particle assemblies is systematically studied using computer simulations of model 
granular packings. Starting from a face-centred cubic array, where all contact forces are identical, an 
increasing number of defects is introduced into the assembly, after which the system is then allowed 
to relax into a new mechanically stable state. Three distinct protocols for imposing disorder are 
compared. A quantitative measure of the disorder is obtained from distributions of the coordination 
number and three-particle contact angle. The distribution of normal contact forces show dramatic 
qualitative changes with increasing disorder. In the regime where the disorder is relatively weak, 
the pressure and the lowest normal mode frequency scale approximately linearly in the coordination 
number, with distance from the crystalline state. These results for P(f) are discussed in the context 
of jamming phenomena in glassy and granular materials. 
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Over the past decade or so it has come to light that the 
way the contact forces are distributed through a static 
pile of grains occurs in a rather inhomogeneous manner 
compared to that of an equilibrium liquid, even though 
the arrangement of the particles are similar in both cases 
. Whereas many pairs of touching particles in a gran- 
ular packing experience only a small force between them, 
others experience much larger forces. In fact, there are 
more of these high-force contacts than one might expect 
if a static packing of particles were truly representative 
of a snapshot of a thermodynamic system. What is be- 
coming increasingly apparent, is that the properties of 
granular packings appear to share many similarities with 
amorphous, glassy phases of their atomic and colloidal 
counterparts Simply due to the fact that grains are 
easier to see than atoms, granular materials now hold a 
prominent role in studies of amorphous systems. 

However, what has yet to be determined in detail is 
how the properties of a static packing depend on the ar- 
rangement of the particles. To date, most studies have 
focused on the two extreme cases of either, fully disor- 
dered particle packings, where there are no long-ranged 
correlations in the particle positions, or perfectly ordered 
arrays. Yet, it is not clear how the properties of a static 
packing change as the configuration is tuned from an or- 
dered array to a disordered state. It is this type of order- 
disorder "transition" that is addressed here. 

With the recent upsurge in the study of granular mate- 
rials, several notable properties of particle packings have 
emerged. One of the most robust, and reproducible, 
features of a disordered grain pile, that sets it apart 
from more traditional crystalline solids, is the probabil- 
ity distribution function of contact forces P(f), where 
/ = FJ < F > is the normal force F between two parti- 
cles in contact, normalized by the configuration-averaged 
force < F >. It has been shown experimentally, numeri- 
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cally, and by computer simulations, that P(f) for a dis- 
ordered particle packing exhibits an exponential'ish Q 
decay at forces /, greater (roughly twice) than the aver- 
age force, and a peak or plateau at small forces (below 
half the average force) 0, 0, IE S • All of these studies 
suggest that P(f) is quite insensitive to system size as 
well dimensionality. Experiments on three dimensional 
packings @ have reported that P(f) is surprisingly insen- 
sitive to the structure of the packing and particle prop- 
erties (such as the friction coefficient). For softer, rubber 
particles there appears to be a cross-over to a power- 
law tail at high forces |10|. Granular dynamics simula- 
tions have shown that for disordered packings, proper- 
ties such as friction and inelasticity, have only a subtle 
effect on P(f) whereas geometrical features of the 
packing, such as coordination number, do depend quite 
sensitively on the parameters chosen Hard-sphere 
simulations of vacancy-diluted crystals show a depletion 
of small forces compared with disordered configurations 

m 

As discussed above, disordered particle packings can 
exhibit a number of subtle differences. However, there is 
one extreme case that is well-defined. For a perfectly or- 
dered, face-centred cubic (fee), crystalline array, all parti- 
cles have the same number of contact neighbours so that 
the coordination number is z = Zf cc = 12, and experience 
the same contact forces, / = 1. The distributions of the 
coordination number P(z), and contact forces P(f), are 
thus known exactly, 

Pfcc(z) = (5(2 - Zf cc ) m 

PfM = S(f-l). ^ 

Clearly, there is a need to extract which properties of a 
static array of particles are responsible for determining 
the way forces are distributed. For this reason, the work 
presented here provides a systematic study as to how 
structural disorder affects the distribution of forces inside 
three dimensional (31?) particle assemblies. 

The computer experiments reported here are for a 
model system: three dimensional, frictionless and non- 
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cohesive, soft-sphere packings, with periodic boundary 
conditions. Two particles are defined to be neighbours 
and interact through a purely-repulsive, short-range, po- 
tential: V(r) — k(d — r) 2 , when their centre-centre sep- 
aration r < d, where d is the sum of their radii, and 
V(r) = 0, otherwise. To create the particle configu- 
rations with varying amounts of disorder, a number of 
protocols were implemented: 

(i) In the first method, N = 16, 384 particles were 
arranged into a 31?, fee array, slightly over- 
compressed to a packing fraction 4> = 0.742, just 
above that of a hard sphere fee array, 0f cc = \/2"k/0>, 
and with z = Z{ cc . Disorder was introduced by ran- 
domly removing, 0.005 < SN < 40% of the parti- 
cles. To remain at constant packing fraction, the 
system was then "quenched" to the same initial 
6 = 0.742, using an energy minimization algorithm 
jl3j. These are denoted quenched packings. The 
resulting changes in the final packings were then 
analysed. In the following figures, each data set for 
a particular SN, was averaged over 5 independent 
realisations. 

(ii) The second protocol was similar to the first, but 
instead of removing particles, particles chosen at 
random had their diameters reduced by 10%, before 
being re-quenched. In this case, the system was 
made increasingly bidisperse. In each of the two 
cases (i) and (ii), the same procedure was again 
repeated using quenched molecular dynamics, with 
N = 16384 and N = 256, 000 particles. 

(iii) The final protocol used started with N — 32, 000 
particles arranged into a fee array at (f){ cc , con- 
strained in the vertical direction by flat base and 
a free, top surface, with periodic boundary con- 
ditions in the horizontal plane [ll|. Disorder was 
then introduced by removing particles at random, 
then allowing the assembly to relax under gravity 
- gravity sedimented packings. 

The various protocols introduced above create configu- 
rations which are, statistically, very similar. Most of the 
results presented here are for the quenched, frictionless 
packings of protocol (i). To clarify the following nomen- 
clature, the configurations generated via protocols (i) and 
(ii) are labelled Cn, where n is an index indicating the 
amount of disorder. Larger n correspond to packings 
with more defects, either by removing particles as in (i), 
or changing the particle size, as in (ii). To provide a 
comparison, two amorphous packings (Al and A2) were 
generated (at two different packing fractions as detailed 
below), using the same interaction potential. 

The average coordination number z, and the three- 
particle contact angle 9, quantify how the protocols gen- 
erate configurations of varying disorder. 6 is defined as 
the angle subtended by particle i and two of its contact 
neighbours: jik. Figure^shows the evolution of the dis- 
tributions P{z). The different configurations are labelled 




FIG. 1: Protocol (i): (a) Coordination number distributions 
P(z), for different amounts of disorder introduced by remov- 
ing a fraction, SN, of the particles. Increasing configuration 
label (C) corresponds to increasing disorder. The "Al" data 
is an amorphous packing generated from an initially random 
distribution of particles quenched to same <f> = 0.742. For 
comparison, the "A2" data is for an amorphous packing at 
random close packing, (j> = 0.64. 

in increasing value according to the number of particles 
removed: configurations C[l-12] correspond to SN/N = 
(0.006, 0.012, 0.018, 0.024, 0.030, 0.06, 0.09, 0.12, 0.15, 
0.18, 0.24, 0.36). The "Al" system was generated from 
an initially random distribution of particles quenched to 
the same 4> = 0.742. The "A2" amorphous configuration 
is close to the random close packing point at </> = 0.64. 

Although the P(z) deviate from the fee delta-function 
even for little disorder, the average coordination z, re- 
main close to Zf cc . The deviation in coordination number 
from the initial crystal state quantitatively characterises 
the disorder. As shown in Fig. [3 for SN << N, 

Sz cx SN, (2) 

where Sz = Z[ cc — z, is the change in the average coor- 
dination number relative to the initial crystal. Thus Sz 
provides a convenient measure of how far away the fi- 
nal configuration is from the crystalline state. Hence, a 
weakly disordered regime can be roughly identified with 
configuratons for which Eq. [21 applies, namely configura- 
tions C[l — 7], for protocol (i). 

Similarly, the distributions in the three-particle con- 
tact angle, P(6), shown in Fig. [31 reflect the fact that 
protocol (i) induces only slight deviations from the fee 
contact topology. Even up to Cll, there are a signifi- 
cant fraction of three-particle collineations (9 = 180°), 
whereas, for amorphous packings these become rare |14| . 

Visualisation of the force networks indicate that the 
forces are very sensitive to the imposed disorder. To illus- 
trate this, Fig. 0] shows the averaged change in the pres- 
sure (normal stresses) between the initial crystal and final 
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FIG. 2: Protocol (i): Relation between the number of defects 
SN, introduced into the packing, and the deviation in the co- 
ordination number 5z, relative to the initial crystal structure. 
A weakly disordered regime can be identified where, Sz oc <5iV. 
Here, this corresponds to configurations C[l — 7]. 
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FIG. 3: Protocol (i): Three-particle contact angle distribu- 
tions for, (top panel) C2 (dashed line) and C6 (full line), and 
(bottom) Clf (dashed) and C12 (full). 



the configurations of two systems with different amounts 
of disorder. 

The distribution of normal contact forces P(f) shown 
in Fig. [S] quantify the resulting changes in the force net- 
work. For a small fraction of defects, upto and including 
configuration C6, the region around the primary peaks 
in P(f) initially undergo Gaussian broadening from the 
original, fee delta-function. Where the peak remains 
identifiable, the standard deviations s, of Gaussian fits 
to the primary P(f) peaks, grow as the location of the 
peaks /peak, move to lower forces with increasing disor- 
der. This is shown in the inset to Fig.l^a). (Recall that 
for the fee array, fp^ ak = 1.) In this regime, the P(f) 
develop additional features either side of / pca k, reflecting 
local force balance in the presence of local defects. These 
features correspond to the regime where Eq. |3 apply. As 



FIG. 4: Normal stress maps of the relative change in pres- 
sure between the initial crystal and the disordered packing. 
Each panel was averaged over 5 realisations. Left: Nearly or- 
dered configuration (C2). Right: intermediate disorder (CI 1). 
Larger stress represented by darker shading. Same greyscale 
in both panels. 




FIG. 5: Protocol (i): Distributions of normal contact forces 
P(f), for varying disorder at fixed cj> = 0.742: (a) linear axes 
and (b) linear- log. In (b), the inner curve is for configuration 
C2 (thick solid line), extending outwards with increasing dis- 
order, C4, C6, C8, C10, and C12 (thick dashed). The circles 
sitting on top of C12 is the data for the "Al" at cj> = 0.742. 
For comparison, the thick dash-dotted line is the "A2" data 
at 4> = 0.64. The inset to (a) shows the standard deviations 
s, of Gaussian fits to the P(f) peaks, growing linearly with 
the peak position /peak, shifted relative to the fee. 



further disorder is introduced, the distributions become 
increasingly asymmetric; showing a dramatic increase in 
the number of very small forces and where the high-/ 
tails become increasingly broad (exponential) |l5j . Thus, 
structure plays a crucial role in determining the hetero- 
geneity of the force network in particle packings. 

For small disorder, protocols (i) and (ii) generate con- 
figurations with properties that are almost indistinguish- 
able. For more disorder, although the configurations gen- 
erated via protocols (i) and (ii) have similar properties, 
there are a few noticeable differences. For this reason, 
the results for protocol (ii) are shown in Fig. |SJ 
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FIG. 6: Protocol (ii): In, (a) the relative change in the co- 
ordination number with respect to the initial monodisperse 
fee array, as the fraction, SN, of particles are made bidis- 
perse. (b) Coordination number distributions P(z), for dif- 
ferent amounts of disorder at fixed (f> = 0.742. Increasing 
configuration label corresponds to increasing disorder. P(6) 
for, (c) C6, and (d) C10 (solid line) and C12 (dashed), (e) 
P(f) over a range of disorder. The inner thick solid line is C2 
and the outer thick line is C12. 



On comparing the results for protocol (i) (Figs. ^-05 
with those of protocol (ii) in Fig. H3 it appears that pro- 
tocol (ii) leads to a more gradual change in the structure, 
compared with protocol (i). For (i), there is a much more 
dramatic change in structure between C10 and C12, than 
for the case of (ii). Still, the general trends remain simi- 
lar, particularly for P(f). 

To provide a comparison with the quenched peri- 
odic packings already presented, results for gravity- 
sedimented, frictionless spheres of protocol (iii), are 
shown in Fig.0 There are a number of subtle differences 
between the coordination number and force distributions 
between protocols (i),(h) and protocol (iii). These differ- 
ences primarily arise as a result of two effects. Firstly, 
the gravity packings are not periodic in all three direc- 
tions, so surface effects play a role. Secondly, the grav- 
ity packings are able to adjust their packing fractions 
with increasing disorder, as shown in Fig. Ufa), whereas 
the fully periodic system are at a fixed packing fraction. 
However, the generic features of the P(z) (Fig.|J{b)) and 
P(f) (Fig. E{c)) distributions are similar. As more dis- 
order is introduced into the lattice structure, the distri- 
butions broaden. In the case of P(f), the high-/ tail 
becomes increasingly exponential. Thus, the fully peri- 
odic systems capture the essential features of the more 
realistic gravity packings. 

The recent focus on P(f) has been emphasised due 
to suggestions that particular properties of P(f) can be 
used to signal the onset of glassiness in a glass-forming 
system, and likewise, the approach of the jammed, static 
state in a granular material or dispersion M ■ 

The con- 
cept being that the development of a peak in P(f), for 
/ < 1, represents the balance of forces required for me- 
chanical stability into the jammed state - development of 
a yield stress. For finite-temperature systems, the jam- 
ming transition of purely repulsive, particles is accompa- 
nied by the development of a peak in P(f) at small /, as 
the temperature is lowered through the glass transition 
temperature |17| . (This picture is not so clear for the 
case of systems with lon ger- range attractive forces, as in 
Lennard- Jones systems [l8L H^. 1 ) At zero-temperature, 
the peak in P(f) flattens into a plateau as the density of 
a jammed packing is lowered towards the jamming tran- 
sition packing fraction 0, Therefore, for purely- 
repulsive, finite-range interactions, this jamming picture 
seems to apply. 

In Fig. [HJ results are shown for two jammed states of 
purely-repulsive systems, each with N = 256000 par- 
ticles. One configuration was generated using a fast 
molecular dynamics quench, from a high (liquid) temper- 
ature to T = 0. The other is a partially melted crystal 
quenched back down to T — 0. Despite the fact that not 
only are both of these systems jammed and that the P(f) 
curves in Fig.[8Ja) sit on top of one another, differences in 
their structures are evident from the radial distribution 
function g(r) and P(9) of Figs. EJb) and (c).This high- 
lights the fact that amorphous and quasi-ordered systems 
can exhibit similar features in the jammed state. 
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FIG. 7: Protocol (iii): Gravity sedimented packings, (a) De- 
pendence of the depth-averaged packing fraction <f>, on SN. 
(b) Coordination number distributions P(z), for different con- 
figurations, labelled with increasing disorder. Configuration 
GA is an amorphous packing, with a packing fraction close 
to the random close packing value « 0.64. (c) and (d) are 
the distributions of normal contact forces P(f), for varying 
disorder on linear axes and linear-log respectively. The thick 
dashed line is GA. 



FIG. 8: (a) The distribution of contact forces P(f) for purely- 
repulsive, soft-spheres at T — and <j> = 0.742, with different 
histories. Soft-sphere glass (solid line) and partially melted 
crystal (dashed line). From P(f) and visualisation of the force 
networks, both configurations appear very similar. Structural 
measures show, however, that the partially melted crystal is 
significantly more ordered than the glassy state: (b) The ra- 
dial distribution function, g(r), shows long-range oscillations 
and, (c) the three-particle contact angle distribution P(6), 
contains additional structure indicative of ordering. 



Increasing disorder affects the mechanical properties of 
the packings subject to external perturbations |2l|. Yet, 
for small amounts of disorder, one expects the configu- 
rations to vary only slightly in their properties from the 



underlying crystal. This is, indeed, the described 
by Eq. [3 To determine how the disorder influences the 
dynamical properties of the packings, the low-frequency 
portion of the "phonon" density of states T>(uS), were ex- 
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FIG. 9: Distribution of the vibrational normal modes of fre- 
quency ui - the density of states - T>(u}), for the quenched 
packings of protocol (i). Only, the Iow-oj region of the dis- 
tribution is shown. With increasing disorder, the number of 
low-frequency modes increases relative to the fee lattice. 



tracted for the N — 16384 systems [23. Figure ED shows 
as a function of frequency u>, for configurations with 
varying amounts of disorder. Increasing disorder leads to 
an increasing population of the low-w region in T>(lo), not 
unlike lattice-disorder models |U • 

The relation, Eq. [3 provides a useful measure of the 
disorder, when the disorder is weak. Figure lTUI shows that 
over this same range in 5z, where the configurations are 
not too different from the original lattice, many of the 
packings' properties scale. In particular, relative to the 
crystal state, the lowest normal mode frequency varies ap- 
proximately linearly with coordination number. This is 
reminiscent of the way jammed amorphous packings be- 
have as the density is lowered towards the random close 
packing point from above H^,|2^|. It is suggestive, there- 
fore, that there may exist a characteristic length scale 
associated with the increasing disorder, though, as yet, 
one has not been identified here. 

In conclusion, the effect of structure on the properties 
of static packings has been studied. Structure plays a 
dramatic role in modifying the distribution of normal 
contact forces P(f) of frictionless particle assemblies, 
from an initial delta-function, for the face-centred cu- 
bic array, to the more familiar 'exponential' decay, with 
finite disorder. Likewise, the distributions of the coor- 
dination number broaden quickly. These two findings 
may, in part, explain the reason why it has been difficult 
to experimentally observe any significant dependence of 
P(f) on structure |9j. Even a relatively small number of 
defects can broaden P(f) quite substantially. 



Dynamical properties of the packings have been inves- 
tigated by extracting statistics on the lowest-lying normal 
modes. As more disorder is imposed, there is an increase 
in the density of states, T>(u)), at small frequency, u). Pro- 
vided the structure remains close to the underlying lat- 
tice, the change in the coordination number relative to 
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FIG. 10: Lowest normal mode frequency Alj (open circles), 
pressure Ap (solid squares), and contact energy Ae (solid cir- 
cles), relative to the original fee values, as a function of the 
degree of disorder, measured by Sz. The solid lines are power- 
laws with exponents one and two. 



the original lattice provides a useful measure of the dis- 
order. When the disorder is weak, vis-a-vis Eq. [3 the 
value of the lowest-lying normal mode frequency scales 
approximately linearly with the disorder. It is some- 
what amusing that these relationships are not too differ- 
ent from what is found in fully disordered packings near 
random close packing |27j . where fractional changes in 
density play a similar role as disorder does here. Between 
these extremes of densely packed ordered arrays and loose 
amorphous packings, lies an intermediate regime that is 
not so well characterised. This intermediate regime re- 
tains a large degree of the contact topology of the origi- 
nal lattice, yet exhibits a strong degree of heterogeneity 
in the contact forces. 
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